/*
 * gmx_meancharges.h
 *
 *  Created on: May 25, 2011
 *      Author: marchi
 */

#ifndef GMX_MEANCHARGES_H_
#define GMX_MEANCHARGES_H_
#include <copyrite.h>
#include <filenm.h>
#include <macros.h>
#include <pbc.h>
#include <smalloc.h>
#include <statutil.h>
#include <vec.h>
#include <xvgr.h>

#include <nbsearch.h>
#include <trajana.h>
#include "gmx_meancharges.h"
#include <fstream>
#include <iostream>
#include <complex>
#include <tpxio.h>
//#include <types/topology.h>
#include "Array.h"
#include "fftw++.h"
#include "rmpbc.h"

using namespace Array;
using namespace fftwpp;

/*! \brief
 * Template analysis data structure.
 */
	typedef struct
	{
		gmx_ana_selection_t *refsel;
		FILE                *fp;
		real                *ave;
		real                *n;
		gmx_ana_nbsearch_t  *nb;
	} t_analysisdata;


	FILE *fp=NULL;
	FILE *fpxplor=NULL;
	t_fileio * topfile;
	real                cutoff = 0;
	gmx_bool            bArg   = FALSE;
	bool bIndex, bCenter=TRUE;
	rvec ngrid={1,1,1};
	int        flags = TRX_READ_X;
	int        flagsb;
    static rvec  trans = {0,0,0};
    static real filter=0.4;
    int pbc_enum;

    enum { ecSel, ecTric, ecRect, ecZero, ecNRR };
    enum
    {
      epSel,
      epNone,
      epComMol,
      epComRes,
      epComAtom,
      epNojump,
      epCluster,
      epWhole,
      epNR
    };

    enum { euSel,euRect, euTric, euCompact, euNR};

	const char *pbc_opt[epNR + 1] =
      { NULL, "none", "mol", "res", "atom", "nojump", "cluster", "whole",
          NULL };
	const char *unitcell_opt[euNR+1] =
      { NULL, "rect", "tric", "compact", NULL };
	const char *center_opt[ecNRR+1] =
      { NULL, "tric", "rect", "zero", NULL };
	int unitcell_enum;
	int ecenter;
	atom_id * cindex, ncent;
	bool bRmPBC,bPBCWhole,bPBCcomRes,bPBCcomMol,bPBCcomAtom,bPBC,bNoJump,bCluster;

	bool bOXT=FALSE,bDipole=FALSE, bXplor=FALSE;
	bool bFilter=FALSE;
	bool bTrans=FALSE;
	const char     *distfile;


static void put_molecule_com_in_box(int unitcell_enum,int ecenter,
                                    t_block *mols,
                                    int natoms,t_atom atom[],
                                    int ePBC,matrix box,rvec x[])
{
    atom_id i,j;
    int     d;
    rvec    com,new_com,shift,box_center;
    real    m;
    double  mtot;
    t_pbc   pbc;

    calc_box_center(ecenter,box,box_center);
    set_pbc(&pbc,ePBC,box);
    if (mols->nr <= 0)
        gmx_fatal(FARGS,"There are no molecule descriptions. I need a tpr file for this pbc option.");
    for(i=0; (i<mols->nr); i++) {
        /* calc COM */
        clear_rvec(com);
        mtot = 0;
        for(j=mols->index[i]; (j<mols->index[i+1] && j<natoms); j++) {
            m = atom[j].m;
            for(d=0; d<DIM; d++)
                com[d] += m*x[j][d];
            mtot += m;
        }
        /* calculate final COM */
        svmul(1.0/mtot, com, com);

        /* check if COM is outside box */
        copy_rvec(com,new_com);
        switch (unitcell_enum) {
        case euRect:
            put_atoms_in_box(box,1,&new_com);
            break;
        case euTric:
            put_atoms_in_triclinic_unitcell(ecenter,box,1,&new_com);
            break;
        case euCompact:
            put_atoms_in_compact_unitcell(ePBC,ecenter,box,1,&new_com);
            break;
        }
        rvec_sub(new_com,com,shift);
        if (norm2(shift) > 0) {
            if (debug)
                fprintf (debug,"\nShifting position of molecule %d "
                         "by %8.3f  %8.3f  %8.3f\n", i+1, PR_VEC(shift));
            for(j=mols->index[i]; (j<mols->index[i+1] && j<natoms); j++) {
                rvec_inc(x[j],shift);
            }
        }
    }
}

static void put_residue_com_in_box(int unitcell_enum,int ecenter,
                                   int natoms,t_atom atom[],
                                   int ePBC,matrix box,rvec x[])
{
    atom_id i, j, res_start, res_end, res_nat;
    int     d, presnr;
    real    m;
    double  mtot;
    rvec    box_center,com,new_com,shift;

    calc_box_center(ecenter,box,box_center);

    presnr = NOTSET;
    res_start = 0;
    clear_rvec(com);
    mtot = 0;
    for(i=0; i<natoms+1; i++) {
#ifdef OLDERVERSION
    	if (i == natoms || (presnr != atom[i].resnr && presnr != NOTSET)) {
#else
        if (i == natoms || (presnr != atom[i].resind && presnr != NOTSET)) {
#endif
        	/* calculate final COM */
            res_end = i;
            res_nat = res_end - res_start;
            svmul(1.0/mtot, com, com);

            /* check if COM is outside box */
            copy_rvec(com,new_com);
            switch (unitcell_enum) {
            case euRect:
                put_atoms_in_box(box,1,&new_com);
                break;
            case euTric:
                put_atoms_in_triclinic_unitcell(ecenter,box,1,&new_com);
                break;
            case euCompact:
                put_atoms_in_compact_unitcell(ePBC,ecenter,box,1,&new_com);
                break;
            }
            rvec_sub(new_com,com,shift);
            if (norm2(shift)) {
                if (debug)
#ifdef OLDERVERSION
                    fprintf (debug,"\nShifting position of residue %d (atoms %u-%u) "
                             "by %g,%g,%g\n", atom[res_start].resnr+1,
                             res_start+1, res_end+1, PR_VEC(shift));
#else
                    fprintf (debug,"\nShifting position of residue %d (atoms %u-%u) "
                             "by %g,%g,%g\n", atom[res_start].resind+1,
                             res_start+1, res_end+1, PR_VEC(shift));
#endif
                for(j=res_start; j<res_end; j++)
                    rvec_inc(x[j],shift);
            }
            clear_rvec(com);
            mtot = 0;

            /* remember start of new residue */
            res_start = i;
        }
        if (i < natoms) {
            /* calc COM */
            m = atom[i].m;
            for(d=0; d<DIM; d++)
                com[d] += m*x[i][d];
            mtot += m;
#ifdef OLDERVERSION
            presnr = atom[i].resnr;
#else
            presnr = atom[i].resind;
#endif
        }
    }
}

static void center_x(int ecenter,rvec x[],matrix box,int n,int nc,atom_id ci[]){
    int i,m,ai;
    rvec cmin,cmax,box_center,dx;

    if (nc > 0) {
        copy_rvec(x[ci[0]],cmin);
        copy_rvec(x[ci[0]],cmax);
        for(i=0; i<nc; i++) {
            ai=ci[i];
            for(m=0; m<DIM; m++) {
                if (x[ai][m] < cmin[m])
                    cmin[m]=x[ai][m];
                else if (x[ai][m] > cmax[m])
                    cmax[m]=x[ai][m];
            }
        }
        calc_box_center(ecenter,box,box_center);
        for(m=0; m<DIM; m++)
            dx[m] = box_center[m]-(cmin[m]+cmax[m])*0.5;

        for(i=0; i<n; i++)
            rvec_inc(x[i],dx);
    }
}

#endif /* GMX_MEANCHARGES_H_ */
